home *** CD-ROM | disk | FTP | other *** search
/ SGI Hot Mix 17 / Hot Mix 17.iso / HM17_SGI / research / lib / map_grid.pro < prev    next >
Text File  |  1997-07-08  |  18KB  |  468 lines
  1. ; $Id: map_grid.pro,v 1.19 1997/03/26 00:15:35 dave Exp $
  2. ;
  3. ; Copyright (c) 1996-1997, Research Systems, Inc.  All rights reserved.
  4. ;    Unauthorized reproduction prohibited.
  5. ;+
  6. ; NAME:
  7. ;    MAP_GRID
  8. ;
  9. ; PURPOSE:
  10. ;       The MAP_GRID procedure draws the graticule of parallels and meridians,
  11. ; according to the specifications established by MAP_SET. MAP_SET must be called
  12. ; before MAP_GRID to establish the projection type, the center of the 
  13. ; projection, polar rotation and geographical limits.
  14. ;
  15. ; CATEGORY:
  16. ;    Mapping.
  17. ;
  18. ; CALLING SEQUENCE:
  19. ;       MAP_GRID
  20. ;
  21. ; INPUTS:
  22. ;    NONE
  23. ;
  24. ; OPTIONAL INPUTS:
  25. ;    NONE
  26. ;    
  27. ; KEYWORD PARAMETERS:
  28. ;
  29. ;   CHARSIZE: The size of the characters used for the labels. The default is 1.
  30. ;      COLOR: The color index for the grid lines.
  31. ;FILL_HORIZON: Fills the current map_horizon.
  32. ;    HORIZON: Draws the current map horizon.
  33. ;  INCREMENT: Determines the spacing between graticle points.
  34. ; GLINESTYLE: If set, the line style used to draw the grid of parallels and 
  35. ;             meridians. See the IDL User's Guide for options. The default is
  36. ;             dotted.
  37. ; GLINETHICK: The thickness of the grid lines. Default is 1.
  38. ;      LABEL: Set this keyword to label the parallels and meridians with their 
  39. ;             corresponding latitudes and longitudes. Setting this keyword to 
  40. ;             an integer will cause every LABEL gridline to be labeled (i.e, 
  41. ;             if LABEL=3 then every third gridline will be labeled). The 
  42. ;             starting point for determining which gridlines are labeled is the
  43. ;             minimum latitude or longitude (-180 to 180), unless the LATS or 
  44. ;             LONS keyword is set to a single value. In this case, the starting
  45. ;             point is the value of LATS or LONS.
  46. ;   LATALIGN: This keyword controls the alignment of the text baseline for 
  47. ;             latitude labels. A value of 0.0 left justifies the label, 
  48. ;             1.0 right justifies it, and 0.5 centers it.
  49. ;     LATDEL: The spacing in degrees between parallels of latitude in the grid. 
  50. ;             If not set, a suitable value is determined from the current map
  51. ;             projection.
  52. ;       LATS: The desired latitudes to be drawn (and optionally labeled). If 
  53. ;             this keyword is omitted, appropriate latitudes will be generated
  54. ;             with consideration of the optional LATDEL keyword. If this
  55. ;             keyword is set to a single value, drawn (and optionally labeled)
  56. ;             latitudes will be automatically generated  as
  57. ;             [...,LATS-LATDEL,LATS,LATS+LATDEL,...] over the extent of the 
  58. ;             map.  The single valued LATS is taken to be the starting point 
  59. ;             for labelling (See the LABEL keyword).
  60. ;     LATLAB: The longitude at which to place latitude labels. The default is 
  61. ;             the center longitude on the map.
  62. ;   LATNAMES: An array specifing the names to be used for the latitude labels.
  63. ;             By default, this array is automatically generated in units of
  64. ;             degrees. This array can be a string array or numeric, but should 
  65. ;             not be of mixed type. When LATNAMES is specified, LATS must also 
  66. ;             be specified, but the number of elments in the two arrays need not
  67. ;             be equal. Extra LATNAMES are ignored, while LATNAMES not supplied
  68. ;             are automatically reported in degrees. LATNAMES can be used when
  69. ;             LATS is set to a single value. It this case, the LATS used will
  70. ;             start at the specified latitude and progress northward, wrapping
  71. ;             around the globe if appropriate. Caution should be used when
  72. ;             using LATNAMES in conjunction with a single LATS, since the
  73. ;             number of visible latitude gridlines is dependent on many factors.
  74. ;   LONALIGN: This keyword controls the alignment of the text baseline for 
  75. ;             longitude labels. A value of 0.0 left justifies the label, 
  76. ;             1.0 right justifies it, and 0.5 centers it.
  77. ;     LONDEL: The spacing in degrees between meridians of longitude in the grid.
  78. ;             If not set, a suitable value is determined from the current map 
  79. ;             projection.
  80. ;     LONLAB: The latitude at which to place longitude labels. The default is 
  81. ;             the center latitude on the map.
  82. ;       LONS: The desired longitudes to be drawn (and optionally labeled). If 
  83. ;             this keyword is omitted, appropriate longitudes will be generated
  84. ;             with consideration of the optional LONDEL keyword. If this 
  85. ;             keyword is set to a single value, drawn (and optionally labeled) 
  86. ;             longitudes will be automatically generated as
  87. ;             [...,LONS-LONDEL,LONS,LONS+LONDEL,...] over the extent of the map.
  88. ;             The single valued LONS is taken to be the starting point for
  89. ;             labeling (See the LABEL keyword).
  90. ;   LONNAMES: An array specifing the names to be used for the longitude labels.
  91. ;             By default, this array is automatically generated in units of
  92. ;             degrees. This array can be a string array or numeric, but should 
  93. ;             not be of mixed type. When LONNAMES is specified, LATS must also 
  94. ;             be specified, but the number of elments in the two arrays need not
  95. ;             be equal. Extra LONNAMES are ignored, while LATNAMES not supplied
  96. ;             are automatically reported in degrees. LONNAMES can be used when
  97. ;             LONS is set to a single value. It this case, the LONS used will
  98. ;             start at the specified longitude and progress eastward, wrapping
  99. ;             around the globe if appropriate. Caution should be used when
  100. ;             using LONNAMES in conjunction with a single LONS, since the number
  101. ;             of visible longitude gridlines is dependent on many factors.
  102. ;ORIENTATION: Specifies the clockwise angle in degrees from horizontal
  103. ;             of the text baseline that the labels should be rotated. Note,
  104. ;             that the rotation used in MAP_SET is also clockwise, but XYOUTS
  105. ;             is normally counterclockwise.
  106. ;OUTPUTS:
  107. ;         Draws gridlines and labels on the current map.
  108. ;
  109. ; EXAMPLE:
  110. ;      map_set,10,20,23.5,/cont,/ortho,/horizon     ; establish a projection
  111. ;         lats=[-65,-52,-35,-20,-2.59,15,27.6,35,45,55,75,89]    
  112. ;                     ; choose the parallels to draw
  113. ;         map_grid,lats=lats,londel=20,lons=-13,label=2,lonlab=-2.5,latlab=7
  114. ;                     ;draw the grid
  115. ;
  116. ; MODIFICATION HISTORY:
  117. ;     Written by:    SVP, May, 23 1996.
  118. ;    DMS, Feb, 1996 Note that this version plots all gridlines
  119. ;            unless a 4 element LIMIT was specified to MAP_SET.
  120. ;       SVP, Nov 1996   Changed over !map1 (new IDL5 maps)
  121. ;    SVP, May 1996    Map_Grid used to live inside of map_set.pro, now
  122. ;                       it lives here.
  123. ;       SVP, May 1996   Changed LABEL to determine the skip between labelled 
  124. ;                       gridlines.
  125. ;                       Also, added the LATS and LONS keywords to give complete
  126. ;                       control over where the labels go.
  127. ;       SVP, Sept 1996  Added LATNAMS,LONAMES, and ORIENTATION keywords
  128. ;    DMS, Nov, 1996    Rev 2 of maps.
  129. ;-
  130.  
  131. function map_grid_incr, span
  132. ;
  133. ; Determine LONDEL or LATDEL if not specified
  134. ;
  135. IF span eq 0 THEN return, 45.
  136. ipow = 0
  137. t = span < 450.
  138. WHILE t lt 5 DO BEGIN t = t * 10 & ipow = ipow +1 & ENDWHILE
  139. increments = [ 1., 2., 4., 5., 10., 15., 30., 45.]
  140. i = 0
  141. WHILE t gt (increments[i] * 10) DO i = i + 1
  142. return, increments[i] / 10^ipow
  143. end
  144.  
  145.  
  146. PRO map_grid, label=label, latdel = latdel, no_grid=no_grid, $
  147.        londel=londel, glinestyle=glinestyle, glinethick=glinethick, $
  148.        lonlab=lonlab, latlab=latlab, lonalign=lonalign, $
  149.        latalign=latalign, lonnames=lonnames, latnames=latnames, $
  150.        whole_map=whole_map, lats=lats, lons=lons, zvalue=zvalue, $
  151.        color=color, t3d=t3d, charsize=charsize, orientation=orientation, $
  152.        horizon=horizon, fill_horizon=fill_horizon, _EXTRA=extra, $
  153.     increment=increment
  154.  
  155.  
  156. ;
  157. ; Put a grid on a previously established map projection.
  158. ;
  159. ; no grid? - in case someone wants just to put labels
  160. no_grid = keyword_set(no_grid)
  161. ;
  162. ; if Label = n, then Labels are added every n gridlines
  163. ;
  164. if n_elements(label) THEN nlabel=fix(abs(label[0])) ELSE nlabel=0
  165. have_lons =  keyword_set(lons)
  166. have_lats =  keyword_set(lats)
  167. have_lonnames = keyword_set(lonnames)
  168. have_latnames =  keyword_set(latnames)
  169. if n_elements(zvalue) eq 0 THEN zvalue = 0
  170.  
  171. ;    Append the graphics keywords:
  172. if n_elements(t3d) then map_struct_append, extra,'T3D',t3d
  173. if n_elements(color) then map_struct_append, extra, 'COLOR',color
  174.  
  175. if n_elements(glinestyle) eq 0 then glinestyle = 1 ;Default = dotted
  176. map_struct_append, extra, 'LINESTYLE', glinestyle ;Append it
  177.  
  178. if n_elements(glinethick) then map_struct_append, extra,'THICK',glinethick
  179. if n_elements(charsize) then map_struct_append, extra,'CHARSIZE',charsize
  180. if n_elements(orientation) then $ ;Orientation is reversed...
  181.   map_struct_append, extra,'ORIENTATION', -1 * orientation
  182.  
  183.  
  184. ;
  185. ; The gridlines can be specified by 
  186. ;
  187. ;  1) an array of lats and/or lons
  188. ;  2) a single lats or lons which is taken to be the center
  189. ;     or 'for sure' lat or lon with gridlines every latdel or londel from it
  190. ;  3) automatically calculated if lats or lons are not specified.
  191. ;
  192. if (!x.type NE 3) THEN $  ; make sure we have mapping coordinates
  193.    message, 'map_grid---Current ploting device must have mapping coordinates'
  194. ;
  195. ; Require that LATS be specified when LATNAMES is ALSO SPECIFIED
  196. ;
  197.  
  198. if have_latnames and not(have_lats) then $
  199.   message,'map_grid---The LATNAMES keyword MUST be used in conjuction '+$
  200.   'with the LATS keyword.'
  201. ;
  202. ; Same for lons
  203. ;
  204. if have_lonnames and not(have_lons) then $
  205.   message,'map_grid---The LONNAMES keyword MUST be used in conjuction '+$
  206.   'with the LONS keyword.'
  207.  
  208. ; Get lat/lon ranges from !MAP. Did MAP_SET specify 4 element limit?
  209. if !map.ll_box(0) ne !map.ll_box(2) or $
  210.   !map.ll_box(1) ne !map.ll_box(3) then begin
  211.     lonmin = !map.ll_box(1) & lonmax = !map.ll_box(3) & lonmax1 = lonmax
  212.     latmin = !map.ll_box(0) & latmax = !map.ll_box(2)
  213. endif else begin                ;If not, use entire globe
  214.     lonmin = -180 & lonmax1 = (lonmax = 180)
  215.     latmin = -90 &  latmax = 90
  216. endelse
  217.  
  218. IF lonmax1 le lonmin THEN lonmax1 = lonmax1 + 360.
  219.  
  220.  
  221.             ;Default grid spacings...
  222. IF n_elements(latdel) eq 0 THEN begin
  223.     latdel = map_grid_incr(latmax-latmin)
  224.     latd = 1
  225. endif else latd = latdel
  226.  
  227. IF n_elements(londel) eq 0 THEN begin
  228.     londel = map_grid_incr(lonmax1-lonmin)
  229.     lond = 1
  230. endif else lond = londel
  231.  
  232.  
  233. ; IF WHOLE_MAP specified, or the deltas are < 1,
  234. ; do not convert the limits into integers
  235. IF (not keyword_set(whole_map)) THEN BEGIN
  236.    IF abs(latmax - latmin) gt 5. and latd ge 1 THEN BEGIN ;Make range integers
  237.     latmin = float(floor(latmin))
  238.     latmax = ceil(latmax)
  239.     ENDIF
  240.    IF abs(lonmax1 - lonmin) gt 5 and lond ge 1 THEN BEGIN
  241.     lonmin = float(floor(lonmin))
  242.     lonmax1 = ceil(lonmax1)
  243.     ENDIF
  244. ENDIF
  245.  
  246. IF N_Elements(label) NE 0 OR N_ELEMENTS(Latlab) ne 0 OR $
  247.   N_Elements(LonLab) NE 0 THEN BEGIN
  248.     IF N_Elements(Latlab) eq 0 THEN Latlab = (lonmin + lonmax1)/2
  249.     IF N_ELements(LonLab) eq 0 THEN LonLab = (latmin +latmax)/2
  250. ENDIF
  251.                                     ; of grid numbers
  252. IF n_elements(latalign) eq 0 THEN latalign = .5    ;Text alignment of lat labels
  253. IF n_elements(lonalign) eq 0 THEN lonalign = .5 ;Text alignment of lon labels
  254. ;
  255. map_proj_info, iproj, CYLINDRICAL=is_cyl, /CURRENT
  256.  
  257. if keyword_set(increment) then step=increment $
  258. else step = 4 < (latmax - latmin)/10.
  259.  
  260. len = long(float((latmax-latmin)) / float(step) + 1.0)
  261. lati = (float(latmax-latmin) / (len-1.)) * findgen(len) + latmin   ;lats
  262.  
  263. ; This fudge avoids curved meridians at the poles because of the split planes
  264. if is_cyl and !map.p0lat eq 0 then begin
  265.     del = 2.0e-2
  266.     if lati(0) eq -90 then lati(0) = del - 90.
  267.     if lati(len-1) eq 90 then lati(len-1) = 90. - del
  268. endif
  269.  
  270.  
  271. First = 1
  272. ;
  273. ; Determine the number of lons and the lon array
  274. ;
  275. CASE n_elements(lons) OF
  276.     0: BEGIN
  277.         n_lons=1+fix((lonmax1-lonmin)/londel)
  278.         longitudes=lonmin+findgen(n_lons)*londel
  279.     ENDCASE
  280.     1: BEGIN
  281.         n_lons=1+fix((lonmax1-lonmin)/londel)
  282.         longitudes=lons[0]+(findgen(n_lons)-n_lons/2)*londel
  283.         index=where(longitudes eq lons[0])
  284.         longitudes=shift(longitudes,-1*index[0])
  285.     ENDCASE
  286.     ELSE: BEGIN
  287.         n_lons=n_elements(lons)
  288.         longitudes=lons
  289.     ENDCASE
  290. ENDCASE
  291.        
  292. ;
  293. ; Build the Longitude Label Flag
  294. ;
  295. CASE nlabel OF
  296.     0: lon_label = bytarr(n_lons)
  297.     ELSE: BEGIN
  298.         lon_label=bytarr(n_lons)
  299.         CASE n_elements(lons) OF
  300.             1: BEGIN
  301.  ; Make sure the center one is set and go out from there
  302.                 index=where(longitudes eq lons[0])
  303.                 i=index[0]      ; Go Up
  304.                 WHILE (i LE n_lons-1) DO BEGIN
  305.                     lon_label[i]=1
  306.                     i=i+nlabel
  307.                 END
  308.                                 ; Go Down
  309.                 i=index[0]-nlabel
  310.                 WHILE (i GE 0) DO BEGIN
  311.                     lon_label[i]=1
  312.                     i=i-nlabel
  313.                 ENDWHILE
  314.             ENDCASE
  315.             ELSE: BEGIN
  316.                                 ; Start with lonmin and label each nlabel point 
  317.                 i=0
  318.                 WHILE (i LE n_lons-1) DO BEGIN
  319.                     lon_label[i]=1 
  320.                     i=i+nlabel
  321.                 END
  322.             END
  323.         ENDCASE
  324.     END
  325. ENDCASE
  326. ;
  327. ;   Dont repeat 180 labelling if the projection is cylindrical and 
  328. ;   Both 180 and -180 are present. This can be defeated by using
  329. ;   LONS=-180
  330. ;
  331. if is_cyl then begin
  332.     id_180 = where(longitudes eq 180,count)
  333.     id_m180 = where(longitudes eq -180,mcount)
  334.     if count gt 0 and mcount gt 0 then begin
  335.         if n_elements(lons) eq 1 then begin
  336.             if lons[0] eq -180 then lon_label[id_180]=0 
  337.         endif else lon_label[id_m180]=0
  338.     endif 
  339. endif
  340.  
  341. ;  Do horizon if specified.
  342. if keyword_set(horizon) then map_horizon, _EXTRA=e
  343. if keyword_set(fill_horizon) then map_horizon, _EXTRA=e, /FILL
  344.  
  345. ;
  346. ; Draw/Label the meridians
  347. ;
  348. FOR i=0,n_lons-1 DO BEGIN
  349.     lon=longitudes[i]
  350.     IF (lon lt -180) THEN lon2 = lon + 360  $
  351.     ELSE IF (lon gt 180) THEN lon2 = -360 +lon $
  352.     ELSE lon2 = lon
  353.     
  354.     IF no_grid eq 0 THEN PLOTS, lon, lati, zvalue, $
  355.       NOCLIP=0, _EXTRA=extra
  356.     
  357.     IF lon2 ne long(lon2) THEN fmt = '(f7.2)' ELSE fmt = '(i4)'
  358.  
  359.     IF lon_label[i] THEN BEGIN
  360.         lonname=strtrim(string(lon2, format=fmt),2)
  361.     ;
  362.     ; replace lonname if a user specified value is available
  363.     ;
  364.     IF have_lonnames then $
  365.      IF i le n_elements(lonnames)-1 then $
  366.       IF (reverse(size(lonnames[i])))[1] eq 7 then $
  367.             lonname=lonnames[i] else $
  368.         lonname=strtrim(string(lonnames[i],format=fmt),2)
  369.     XYOUTS, lon, LonLab, ALIGNMENT=lonalign,lonname, NOCLIP=0, $
  370.           Z=zvalue, _EXTRA=extra 
  371.     ENDIF
  372. ENDFOR
  373.  
  374.     step = 4 < (lonmax1 - lonmin)/10.
  375. len = (lonmax1-lonmin)/step + 1
  376. loni = findgen(len)*step + lonmin
  377.  
  378. IF (loni[len-1] NE lonmax1) THEN  BEGIN
  379.     loni = [loni, lonmax1]
  380.     len = len + 1
  381. ENDIF
  382.  
  383. ;
  384. ; Determine the number of lats and the lat array
  385. ;
  386. CASE n_elements(lats) OF
  387.         0: BEGIN
  388.             n_lats=1+fix((latmax-latmin)/latdel)
  389.             latitudes=latmin+findgen(n_lats)*latdel
  390.            END
  391.     1: BEGIN
  392.             n_lats=1+fix((latmax-latmin)/latdel)
  393.             latitudes=lats[0]+(findgen(n_lats)-n_lats/2)*latdel
  394.             index=where(latitudes eq lats[0])
  395.             latitudes=shift(latitudes,-1*index[0])
  396.            END
  397.      ELSE: BEGIN
  398.              n_lats=n_elements(lats)
  399.              latitudes=lats
  400.            END
  401. ENDCASE
  402.  
  403. ;
  404. ; Build the Longitude Label Flag
  405. ;
  406. CASE nlabel OF
  407.         0: lat_label=bytarr(n_lats) > 0
  408.      ELSE: BEGIN
  409.            lat_label=bytarr(n_lats) > 0
  410.             CASE n_elements(lats) OF
  411.                  1: BEGIN
  412.             ; Make sure the center one is set
  413.                     ; and go out from there
  414.             ;
  415.                        index=where(latitudes eq lats[0],count)
  416.                     ; Go Up
  417.                        i=index[0]
  418.                        WHILE (i LE n_lats-1) DO BEGIN
  419.                          lat_label[i]=1
  420.                          i=i+nlabel
  421.                        END
  422.                     ; Go Down
  423.                        i=index[0]-nlabel
  424.                        WHILE (i GE 0) DO BEGIN
  425.                          lat_label[i]=1
  426.                          i=i-nlabel
  427.                        END
  428.                     END
  429.               ELSE: BEGIN
  430.             ; Start with latmin and label each nlabel point 
  431.                     ; 
  432.                         i=0
  433.                         WHILE (i LE n_lats-1) DO BEGIN
  434.                              lat_label[i]=1 
  435.                  i=i+nlabel
  436.                         END
  437.                      END
  438.             
  439.         ENDCASE
  440.            END
  441. ENDCASE
  442.  
  443. ;
  444. ; Draw/Label the parallels of latitude
  445. ;
  446. FOR i=0,n_lats-1 DO BEGIN
  447.    lat=latitudes[i]
  448.    IF lat ne long(lat) THEN fmt = '(f7.2)' ELSE fmt = '(i4)'
  449.    IF lat_label[i] THEN BEGIN
  450.       latname=strtrim(string(lat, format=fmt),2)
  451.       ; 
  452.       ; replace latname if a user specified value is available
  453.       ;
  454.       IF have_latnames then $
  455.        IF i le n_elements(latnames)-1 then $
  456.         IF (reverse(size(latnames[i])))[1] eq 7 then $
  457.         latname=latnames[i] else $
  458.             latname=strtrim(string(latnames[i],format=fmt),2)
  459.         XYOUTS, latlab, lat, latname, alignment=latalign, NOCLIP=0, $
  460.           Z=zvalue, _EXTRA=extra
  461.    ENDIF
  462.   
  463.    IF (no_grid eq 0) and (abs(lat) ne 90) THEN PLOTS, loni, lat, zvalue, $
  464.     NOCLIP=0, _EXTRA=extra
  465. ENDFOR
  466.  
  467. end
  468.